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Abstract 

English version: We are interested in simulating blood flow in arter- 
ies with variable elasticity with a one dimensional model. We present a 
well-balanced finite volume scheme based on the recent developments in 
shallow water equations context. We thus get a mass conservative scheme 
which also preserves equilibria of Q = 0. This numerical method is tested 
on analytical tests. 

Version Frangaise : Nous nous interessons a la simulation d'ecoulements 
sanguins dans des arteres dont les parois sont a elasticity variable. Ceci 
est modelise a l'aide d'un modele unidimensionnel. Nous presentons un 
schema "volume fini equilibre" base sur les developpements recents ef- 
fectues pour la resolution du systeme de Saint- Venant. Ainsi, nous obtenons 
un schema qui preserve le volume de fluide ainsi que les equilibres au repos: 
Q = 0. Le schema introduit est teste sur des solutions analytiques. 



Introduction 

We consider the following system of mass and momentum conservation with non 
dimensionless parameters and variables, which is the ID model of blood flow in 
an artery or a vessel with non uniform elasticity (it is rewritten in a conservative 
form compared to what we usually find in litterature) 



d t A + d x Q = 
d t Q + d x 



kA 3/2 



^ (d x Ao - ^VAd x k 



Q , (i) 



.4 



with Aq = k^/Ao and where A(x, t) is the cross-section area (A = irR 2 with R 
the radius of the arteria), Q(x, t) = A(x, t)u(x, t) the flow rate or the discharge, 
u(t,x) the mean flow velocity p the blood density, Aq(x) the cross section at 
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rest and k(x) the stiffness of the artery. System (1) is into the form of the Saint- 
Venant problem with variable pressure presented in [3] . We have to mention that 
arterial pulse wavelengths are long enough to justify the use of a ID model rather 
than a 3D model when a global simulation of blood flow in the cardiovascular 
system is needed. 



1 Numerical method 

Since [2, 8], it is well known (in the shallow water community) that the scheme 
should be well-balanced for good source term treatment, i.e. the scheme should 
preserve at least some steady states. For system (1), we should preserve at least 
the " man at eternal rest" or " dead man equilibrium" [6] (without artifacts such 
as [10]), it writes 

u = 

l A V2 k _A Ao = Cst , (2) 

this means that steady states at rest are preserved (this is the analogous of the 
"lake at rest" equilibrium). Thus we use the scheme proposed in [3, p. 93-94] for 
that kind of model. This is a finite volume scheme with a modification of the 
hydrostatic reconstruction (introduced in [1, 3] for the shallow water model). 

1.1 Convective step 

For the homogeneous system 

dtU + d x F(U, Z) = 0, (3) 

which is (1) with: 

^=(g)' z =(t) Md F & k) = ( QVA + J'vwcp) 



an explicit first order conservative scheme writes 

U i - U t , r i+l/2 r i-l/2 

At Ax 
where [/" is an approximation of U 



= 0, (4) 



i r x i+i/2 
f/f ~ — / U(x,t n )dx. 



'Xt-1/2 

i refers to the cell C 2 ; = (ajj-i^jiEt+i/a) = (#1-1/21 ^i- 1/2 + Ax) and n to time 

t„ with t n+ i - t„ = At. 

The two points numerical flux 

= UjXi, ^+1/2)1 

with k* +1 i 2 = max(fc, , fc J+ i), is an approximation of the flux function F(U,Z) 
at the cell interface i + 1/2. This numerical flux will be detained in subsection 
1.3. 
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1.2 Source terms treatment 

In system (1), the terms A(d x Ao — 2\ /r Ad x k/3)/(y/np) are involved in steady 
states preservation, they need a well-balanced treatment: the variables are re- 
constructed locally thanks to a variant of the hydrostatic reconstruction [3, 
p.93-94] 

y / A+i/2L = ma,x(k iy ^A~ + min(A.4oi+i/2, 0), 0)/fc* +1/2 

Ui+i/2L = {Aj + i/2LiA i+ ij2L-U i) t . . 

\/A+i/2R = max.(k i+ly /A i + 1 - max(A^loi+i/2, 0), 0)/fc* +1/2 ' 

, Ui+\/2R = {A i+ i/2R, A i+ i/2R-Ui + iY 



with A^loi+1/2 = Aoi+i-Aoi = k i+ i y/A Qx+1 - k t y/A~^~ and k* +1/2 = max(fc 4 , k i+1 ). 
For consistency, the scheme (4) is modified as follows 

= U i~^ { F i + 1/2L - F? + 1/2 R ) , (6) 



rpn 77171 j o 

P i+1/2L ~ ^i+1/2 + a i+l/2L 

F i-1/2R = ^i-1/2 + S t-1/2R 

F i+l/2 = F (Ut+1/2L, U l+1 /2R, k* +1/2 




where 



with 



+i/2L-y ■p { A?,k t )-V(A? +1/2L ,k; +1/2 ) 

Si -^ R= ( P(A?,h)-V(%_ 1/2R ,k*_ 1/2 ) 

and V(A, k) = kA 3 ^ 2 / (3py/n). Thus the variation of the radius and the varying 
elasticity are treated under a well-balanced way. In system (1), the friction term 
—CfQ/A is treated semi-implicitly. This treatment is classical in shallow water 
simulations [4, 11] and had shown to be efficient in blood flow simulation as well 
[G]. This treatment does not break the "dead man" equilibrium. It consists in 
using first (6) as a prediction step without friction, i.e.: 

U i =U ? ( F i+i/2£ _ F T-i/2r) > 

then we apply a semi- implicit friction correction on the predicted values (U*): 

Thus we get the corrected velocity it™ +1 and we have A™ +1 = A*. 
1.3 HLL numerical flux 

As presented in [6], several numerical fluxes might be used (Rusanov, HLL, 
VFRoe-ncv and kinetic fluxes). In this work we will use the HLL flux (Harten 
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Figure 1: No spurious flows (right) are generated by a change of elasticity (left). 



Lax and van Leer [9] ) because it is the best compromise between accuracy and 
CPU time consuming (see [5, chapter 2]). It writes: 

F(U L ,k*) if0<ci 

c 2 F(U L ,k*) - ciF(U R ,k*) cxc-2 , rT \ - f 

T{\J h ,U R ,k ) = < 1 (U R -U L ) ifci<0<c 2 , 

c 2 - ci c 2 - ci 

F(U R ,k*) ifc 2 < 

with 

ci = inf ( inf X 3 ;(U, k*)) and c 2 = sup ( sup Xj(U,k*)), 
u=u L ,u R je{i,2} u=u L ,u R j<e{i : 2} 

where Xi(U,k*) and A 2 (C7, k*) arc the eigenvalues of the system and k* — 
max(fc£, k R ). 

To prevent blow up of the numerical values, we impose the following CFL 
(Courant, Friedrichs, Levy) condition 

Ax 

At < n C FL 



max(|ui| + Ci) ' 



where ci = \/ kiy/Ai '/ (2p^/n) and ticfl = 1- 



2 Some numerical results 

2.1 "The stented man at eternal rest" 

In this test, we consider a configuration with no flow and with a change of artery 
elasticity k(x), this is the case for a dead man with a stented artery (see Figure 
1 left). The section of the artery is constant Rq(x) = 4.0 10 -3 m and the velocity 
is u(x 7 t) = Om/s. We use the following numerical values: J = 50 cells, Cf = 0, 
p = 1060m 3 , L = 0.14m, T en d = 5s. As initial conditions, we take a fluid at rest 
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Q(x,0) = 0m 3 /s and 
fc 



Ak . 
fc + — — sin 



fc + Ak 
Ak 



cos 



X — X\ 

x 2 - x\ 

x — X3 
Xi - x 3 ' 



7T - n/2 + 1 



1 



if x G [0 : xi] U [x 4 ■■ L] 
if x g]xx : x 2 [ 
if x G [x 2 : x 3 ] 
if x G]x3 : Xi[ 



with fc = 1.0 10 8 Pa/m, Afc = 6.0 10 7 Pa/m, Xi = 1.0 10~ 2 m, x 2 = 3.05 10~ 2 m, 
x 3 = 4.95 10~ 2 m and x 4 = 7.0 10" 2 m. 

The steady state at rest is perfectly preserved in time, we do not notice any 
spurious oscillation (see Figure 1 right). 



2.2 Wave reflection-transmission in a stented artery 



We now observe the reflexion and transmission of a pulse through a sudden 
change of artery elasticity (from kfi to fci with fc^ > /cr) in an elastic tube of 
constant radius (sec Figure 2 left). We take the following numerical values: J = 
1500 cells, C f = 0, k L = 1.6 10 8 Pa/m, k R = 1. 10 8 Pa/m, Afc = 6. 10 7 Pa/m, p = 
1060m 3 , Rq = 4.0 10" 3 m, L = 0.16m, T end = 8.0 10~ 3 s, c L = v / k L R/{2p) ~ 
17.37m/s and cr = y/ k^R/(2p) ~ 13.74m/s. We take a decreasing elasticity 
on a rather small scale: 




cos 



x 2 - xi 



if x G [0 : xx] 
if x G]xx : X2] 
else 



with xx = 19L/40 and x 2 = L/2. As initial conditions, we consider a fluid at 
rest Q(x,0) = 0m 3 /s and the following perturbation of radius: 



i?(x,0) 



1 



/100 / 65L 
3m — — 7T x — 



\20L 



100 



if x G [65L/100 : 85L/100] 
else 



Ro(x) 
Ro{x) 

with e = 1.0 10 -2 . The expression for the pressure is 

p{x, t) = p Q (x) + k(x)(R(x, t) - Ro(x)), 

where po is the external pressure. 

The numerical results perfectly match with the predictions for a linearized 
flow. We get the predicted amplitudes both for the transmitted and the reflected 
waves (see Figure 2 right). 



2.3 Wave "damping" 

In this case, the elasticity is constant in space. We consider the viscous term in 
the linearized momentum equation. A periodic signal is imposed at inflow as a 
perturbation of a steady state (Rq = Cst, uq = 0) with a constant section at 
rest. We take R = Rq+eRi and u = 0+eux, where (Ri,m) is the perturbation of 
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Figure 2: Across a discontinuity of k (left), an initial pulse evolves (right) ; it 
is transmitted and reflected. 



the steady state. Looking for progressive waves {i.e. under the form e^ wt Kx '), 
we take for the imposed incoming discharge 

Qb(t) = Q(t,x = 0) = Qamp sin(wt) m 3 /s. 

Thus, we have a damping wave in the domain 

f if k r x > cut 

\ Qamp sin(u;i — k r x)e kiX if k r x < tut ' 

with ui = 2n /T pu i S e, T pu i se the time length of a pulse and K = k r + iki the wave 
vector. 

We use the following numerical values: J = 750 cells, p = 1060m 3 , L = 3m, 
k = 1.10 8 Pa/m, i? = 4.10~ 3 m and T end = 5s. We consider both C/ = and 
Cf = 0.005053. As initial conditions, we take a fluid at rest Q(x, 0) = 0m 3 /s 
and as input boundary condition 

Qb(t) = Q amp sm(ujt), 

with uj = 2Tr/T pu i se = 27r/0.5s and Q am p = 3.45 10~ 7 m 3 /s. The output is an 
outgoing wave. 

The results are closed to the analytical solution (see Figure 3). We notice a 
small numerical diffusion for C/ = 0.005053. 



3 Conclusion 

In this work, we have considered the ID model of flow in an artery with varying 
elasticity and constant section. We have presented a well-balanced finite volume 
scheme. Thus we get a mass conservative scheme. Moreover, the well-balanced 
property allows to have a good treatment of the source, i.e. we do not get 
numerical artifacts. This numerical method gave good results on numerical 
tests. In future works, we will have to add some extra source terms in order 
to get a more realistic model. These extra terms will require to develop a low 
diffusive high order scheme in the spirit of [7]. Moreover, this will improve the 
accuracy of the scheme. And we will also have to test more complex cases such 
as bifurcations and networks. 
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